******************************************************************************
*										
*                           Boese & Eberhardt							
*            
*                          Facets of Democracy
*													
******************************************************************************
*                   PCDID - heterogeneous model (1949-2018)
*			    (Specification with exports/trade & pop growth)
* 				  Mid-level indices: Interaction Effects
*            Using the alternative model (f^A, f^B instead of f^AB)
******************************************************************************
*                       Using standardized thresholds 
******************************************************************************
*                           Created: 17th July 2024					
*                      Last Changed: 18th July 2024					
******************************************************************************
*                              Markus Eberhardt						
******************************************************************************
*                 School of Economics, University of Nottingham,			
*                  University Park, Nottingham NG7 2RD, England			
*                     email: markus.eberhardt@nottingham.ac.uk				
*                   web: https://sites.google.com/site/medevecon/			
******************************************************************************

clear all
set matsize 11000

***
*** Paths
***

global path "D:/Dropbox/Facets of Democracy"
global rawpath "D:/Dropbox"
global path "/Users/lezme/Dropbox/Facets of Democracy"
global rawpath "/Users/lezme/Dropbox"

global data "$path/Data"
global dofolder "$path/Do_files"
global otherdata "$path/Data"
global outputfolder "$path/Output"
global texfolder "$rawpath/Apps/Overleaf/Facets of Democracy"

global xtmgpath = "D:/Dropbox/Literature/econometric issues/Nonstationary Panel Metrics/Stata"
global xtmgpath = "/Users/lezme/Dropbox/Literature/econometric issues/Nonstationary Panel Metrics/Stata"

****
**** Get merged dataset and prepare variables
****

use "$data/madd_ert_vdem_2021.dta", clear
tsset nwbcode year

* Maddison: per capita GDP (in logs)*100 and population growth rate
gen y = ln(gdppc)*100
gen lpop = ln(pop)
gen dlpop = 100*d.lpop

* DOTS data: export/trade
gen ex = ex_trade*100

xtreg y dlpop ex v2x_libdem if year>=1959, fe

drop csum 
gen c=1 if e(sample)
by nwbcode: egen csum=sum(c) if c==1
gen sample= (e(sample))

* V-Dem: democracy dummies from high and mid-level indicators
local z "libdem polyarchy liberal"
foreach k of local z{
	sum v2x_`k' if sample==1
	scalar mean`k' = r(mean) 
	scalar sd`k' = r(sd) 
	gen st_`k' = (v2x_`k'-mean`k')/sd`k' if sample==1
}

* Positive values
local z "0 25 50"
foreach k of local z{
	local norm_`k' = `k'/100
	gen dem_ld_`k' = (st_libdem>=`norm_`k'') & !missing(v2x_libdem)
	label var dem_ld_`k' "Libdem Threshold > mean + .`k' sd"
	gen dem_poly_`k' = (st_polyarchy>=`norm_`k'') & !missing(v2x_polyarchy)
	label var dem_poly_`k' "Polyarchy Threshold > mean + .`k' sd"
	gen dem_lib_`k' = (st_liberal>=`norm_`k'') & !missing(v2x_liberal) 
	label var dem_lib_`k' "Liberal Component Threshold > mean + .`k' sd"
}

* Negative values
local z "50 25"
foreach k of local z{
	local norm_neg`k' = -`k'/100
	gen dem_ld_neg`k' = (st_libdem>=`norm_neg`k'') & !missing(st_libdem)
	label var dem_ld_neg`k' "Libdem Threshold > mean - .`k' sd"
	gen dem_poly_neg`k' = (st_polyarchy>=`norm_neg`k'') & !missing(st_polyarchy)
	label var dem_poly_neg`k' "Polyarchy Threshold > mean - .`k' sd"
	gen dem_lib_neg`k' = (st_liberal>=`norm_neg`k'') & !missing(st_liberal) 
	label var dem_lib_neg`k' "Liberal Component Threshold > mean - .`k' sd"
}


* Positive values
local z "125"
foreach k of local z{
	local norm_`k' = `k'/1000
	gen dem_ld_`k' = (st_libdem>=`norm_`k'') & !missing(v2x_libdem)
	label var dem_ld_`k' "Libdem Threshold > mean + .`k' sd"
	gen dem_poly_`k' = (st_polyarchy>=`norm_`k'') & !missing(v2x_polyarchy)
	label var dem_poly_`k' "Polyarchy Threshold > mean + .`k' sd"
	gen dem_lib_`k' = (st_liberal>=`norm_`k'') & !missing(v2x_liberal) 
	label var dem_lib_`k' "Liberal Component Threshold > mean + .`k' sd"
}


* Negative values
local z "125"
foreach k of local z{
	local norm_neg`k' = -`k'/1000
	gen dem_ld_neg`k' = (st_libdem>=`norm_neg`k'') & !missing(st_libdem)
	label var dem_ld_neg`k' "Libdem Threshold > mean - .`k' sd"
	gen dem_poly_neg`k' = (st_polyarchy>=`norm_neg`k'') & !missing(st_polyarchy)
	label var dem_poly_neg`k' "Polyarchy Threshold > mean - .`k' sd"
	gen dem_lib_neg`k' = (st_liberal>=`norm_neg`k'') & !missing(st_liberal) 
	label var dem_lib_neg`k' "Liberal Component Threshold > mean - .`k' sd"
}


save "$data/madd_ert_vdem_2021_helper_inter.dta", replace
 
use "$data/madd_ert_vdem_2021_helper_inter.dta", clear
 
gen dem_inter_0 = dem_poly_0 * dem_lib_0 if sample==1
tab  dem_poly_0 dem_lib_0 if sample==1, cell
* 10.8%
tab  dem_inter_0 dem_poly_0 if sample==1, cell
* 3.2%
tab  dem_inter_0 dem_lib_0 if sample==1, cell
* 7.6% 

xtlogit dem_inter_0 y if sample==1, fe iterate(1)
tab  dem_poly_0 dem_lib_0 if e(sample), cell

reg dem_poly_0 i.nwbcode dem_inter_0 dem_lib_0 if e(sample)
drop dem_inter_0

****
**** PCDID 
**** 

estimates clear 
cls 

local i = 1

**** Select the variables, sample period, factor-augmentation 
global xvars1 		"ex dlpop" // X variables used to construct factors 
global xvars2 		"ex dlpop" // Xs included in the PCDID (leave blank for 'plain vanilla' model)
global factors 		"6"		 // max # of factors (over-prediction is not as bad as underprediction)
global uhat1 		"uhatMG1" // uhatMG: residuals from heterog regression (first control group)
global uhat2 		"uhatMG2" // uhatMG: residuals from heterog regression (second control group)
global uhatint 		"uhatMGint" // uhatMG: residuals from heterog regression (interaction control group)

*** Pick which threshold definition is used (mean or mean +- sd)
*local sdvari "neg25 neg125 0 125 25"
local sdvari "neg25 neg125 0 125 25"
* 0 is the model using the standardised mean of the democracy index
* the other indices are neg25 neg125 125 25
foreach ksdvari of local sdvari{
		
	*** List democracy indicators employed here (list1 always poly or lib, list_1 the other one)
	global demo_list1	"lib_`ksdvari'" 
	global demo_list2	"dem_lib_`ksdvari'" 
	global demo_list_1	"poly_`ksdvari' " 
	global demo_list_2	"dem_poly_`ksdvari'" 

	run "$xtmgpath/xtmg.ado"

	* Start of loop for different end years 
	local zzz "2018"
	foreach kkk of local zzz{
		global endyear 		"`kkk'"		// 2018

		* Start of loop for different start years 
		local zz "1959"
		foreach kk of local zz{
			global startyear 	"`kk'"	// 1959

			* Start of loop for the first democracy indicators (e.g. poly) 
			local z "$demo_list1"
			foreach k of local z{
				global democracy 	"`k'"	// picks one democracy indicator

				* Start of loop for the second democracy indicator (e.g. lib) !!! Underscore_1 !!!
				local z2 "$demo_list_1"
				foreach k2 of local z2{
					global democracy_too 	"`k2'"	// picks second democracy indicator (from a list)

					* Start of factor loop: how many factors for dem_list1 variable to be included in the PCDID
					forvalues f=1/$factors{
					preserve
					quietly{
						
						tsset nwbcode year
						keep if year>=$startyear & year<=$endyear
						
						gen dem_dum = dem_$democracy
						gen dem_dumm = dem_$democracy_too
						gen dem_inter = dem_dum * dem_dumm
						sort nwbcode year
						
						* Cross the threshold
						gen count_$democracy = 1 if ((dem_dum==0 & l.dem_dum==1) | (dem_dum==1 & l.dem_dum==0))
						gen count_$democracy_too = 1 if ((dem_dumm==0 & l.dem_dumm==1) | (dem_dumm==1 & l.dem_dumm==0))
						gen count_inter = count_$democracy * count_$democracy_too
						
						by nwbcode: egen dem_dum_count_$democracy = sum(count_$democracy)
						by nwbcode: egen dem_dumm_count_$democracy_too = sum(count_$democracy_too)
						by nwbcode: egen dem_dum_count_inter = sum(count_inter)
						
						replace count_$democracy = 0 if (dem_dum==0 & l.dem_dum==1)
						replace count_$democracy_too = 0 if (dem_dumm==0 & l.dem_dumm==1)
						replace count_inter = 0 if ((dem_dumm==0 & l.dem_dumm==1) | dem_dum==0 & l.dem_dum==1)
						
						* Cross the threshold from below
						by nwbcode: egen trudem_dum_count_$democracy = sum(count_$democracy)
						by nwbcode: egen trudem_dumm_count_$democracy_too = sum(count_$democracy_too)
						by nwbcode: egen trudem_dum_count_inter = sum(count_inter)

						by nwbcode: egen dem_dum_exposure_$democracy = sum(dem_dum)
						by nwbcode: egen dem_dumm_exposure_$democracy_too = sum(dem_dumm)
						by nwbcode: egen dem_dum_exposure_inter = sum(dem_inter)
	
						by nwbcode: egen sample_obs_$democracy = sum(sample)
						by nwbcode: egen sample_obs_$democracy_too = sum(sample)
						by nwbcode: egen sample_obs_inter = sum(sample)

						gen autoc_exposure_$democracy = sample_obs_$democracy - dem_dum_exposure_$democracy
						gen autoc_exposure_$democracy_too = sample_obs_$democracy_too - dem_dumm_exposure_$democracy_too
						gen autoc_exposure_inter = sample_obs_inter - dem_dum_exposure_inter

						*** 1 ***
						* Sample identifiers for the first dem variable
						xtlogit dem_dum y if sample==1, fe iterate(1)
						* Treat are those which experienced crossing the threshold
						gen treat = (e(sample))
						replace treat = . if sample==0
						* Below are those which are always below the threshold
						gen below = (!e(sample) & sample==1 & !missing(dem_$democracy))
						replace below = . if sample==0
						* Always above threshold (adjust)
						gen above = (below==1 & dem_$democracy ==1)
						replace above = . if sample==0
						* Always below threshold (adjust)
						replace below = . if below==1 & dem_$democracy==1
						
						*** 2 ***
						* Sample identifiers for the second dem variable
						xtlogit dem_dumm y if sample==1, fe iterate(1)
						* Treat are those which experienced crossing the threshold
						gen treat_too = (e(sample))
						replace treat_too = . if sample==0
						* Below are those which are always below the threshold
						gen below_too = (!e(sample) & sample==1 & !missing(dem_$democracy_too))
						replace below_too = . if sample==0
						* Always above threshold (adjust)
						gen above_too = (below_too==1 & dem_$democracy_too ==1)
						replace above_too = . if sample==0
						* Always below threshold (adjust)
						replace below_too = . if below_too==1 & dem_$democracy_too==1

						*** 3 ***
						* Sample identifiers for the interaction variable
						gen treat_inter = treat_too*treat
						gen below_inter = below_too*below
						gen above_inter = above_too*above
						
						*** 1 ***
						*** Estimate linear heterog regression of y on x AND second dem dummy in the first control group and collect residuals
						xtmg y dem_dumm $xvars1 if below==1 & sample==1 & above==0, res(uhatMG1)
						local below = e(N_g)
						local belowobs = e(N)
						*** Estimate linear pooled regression of y on x in the control group and collect residuals
						xtreg y dem_dumm $xvars1 if below==1 & sample==1 & above==0, fe
						predict uhatP1, e
						*** Compute residual period average in control group 
						sort year 
						by year: egen uCT1 = mean($uhat1)
						sort nwbcode year							
						*** Extract factors from the residuals - decide which version ($uhat) to use
						xtreg $uhat1 if sample==1, fe
						gen t_max = e(g_max)
						gen N_max = e(N_g)
						if t_max > N_max{
							regife $uhat1 if sample==1, f(cfactor = year floading = nwbcode, `f') maxiterations(1)
						}
						else if  N_max > t_max {
							regife $uhat1 if sample==1, f(floading = nwbcode cfactor = year, `f') maxiterations(1)
						}
						drop t_max N_max
						*** Create factors for all countries
						forvalues l=1/`f'{
							sort year
							by year: egen c1factor`l'N = mean(cfactor`l')		
							drop cfactor`l' 
						}
						drop floading*
						* Control group 1: c1factor*N
						
						*** 2 ***
						*** Estimate linear heterog regression of y on x AND first dem dummy in the second control group and collect residuals
						xtmg y dem_dum $xvars1 if below_too==1 & sample==1 & above_too==0, res(uhatMG2)
						local below_too = e(N_g)
						local belowobs_too = e(N)
						*** Estimate linear pooled regression of y on x in the control group and collect residuals
						xtreg y dem_dum $xvars1 if below_too==1 & sample==1  & above_too==0, fe
						predict uhatP2, e
						*** Compute residual period average in control group 
						sort year 
						by year: egen uCT2 = mean($uhat2)
						sort nwbcode year					
						*** Extract factors from the residuals - decide which version ($uhat) to use
						xtreg $uhat2 if sample==1, fe
						gen t_max = e(g_max)
						gen N_max = e(N_g)
						if t_max > N_max{
							regife $uhat2 if sample==1, f(cfactor = year floading = nwbcode, `f') maxiterations(1)
						}
						else if  N_max > t_max {
							regife $uhat2 if sample==1, f(floading = nwbcode cfactor = year, `f') maxiterations(1)
						}
						drop t_max N_max
						*** Create factors for all countries
						forvalues l=1/`f'{
							sort year
							by year: egen c2factor`l'N = mean(cfactor`l')		
							drop cfactor`l' 
						}
						drop floading*
						* Control group 2: c2factor*N
						
						* NOT NECESSARY IN THIS SPECIFICATION
						*** 3 ***
						*** Estimate linear heterog regression of y on x in the interaction control group and collect residuals
						xtmg y $xvars1 if below_inter==1 & sample==1 & above_inter==0, res(uhatMGint)
						local below_inter = e(N_g)
						local belowobs_inter = e(N)
						*** Estimate linear pooled regression of y on x in the control group and collect residuals
						xtreg y $xvars1 if below_inter==1 & sample==1  & above_inter==0, fe
						predict uhatPint, e
						*** Extract factors from the residuals - decide which version ($uhat) to use
						xtreg $uhatint if sample==1, fe
						gen t_max = e(g_max)
						gen N_max = e(N_g)
						if t_max > N_max{
							regife $uhatint if sample==1, f(cfactor = year floading = nwbcode, `f') maxiterations(1)
						}
						else if  N_max > t_max {
							regife $uhatint if sample==1, f(floading = nwbcode cfactor = year, `f') maxiterations(1)
						}
						drop t_max N_max
						*** Create factors for all countries
						forvalues l=1/`f'{
							sort year
							by year: egen c3factor`l'N = mean(cfactor`l')		
							drop cfactor`l' 
						}
						drop floading*
						* Control group 3 for interaction: c3factor*N
						
						gen dem_interaction = dem_$democracy * dem_$democracy_too
						
						*** Test of controls ***
						xtlogit dem_interaction y $xvars2 ///
							if sample==1 & treat==1, fe iterate(1)
						xtmg dem_interaction $xvars2 c1factor*N c2factor*N ///
							if e(sample)
						testparm $xvars2
						local chi2_stat = r(chi2)
						local chi2_p = r(p)

						*** Alpha test auxiliary regressions ***
						xtmg y uCT1 uCT2 ///
							dem_interaction ///
							$xvars2 ///
							if sample==1 & treat_inter==1
						mat alphas = e(betas)
						svmat alphas, name(alphas)
						gen alpha01 = alphas1
						replace alpha01 = . if (alphas3)==0
						mkmat alpha01, mat(alpha01)
						gen alpha02 = alphas2
						replace alpha02 = . if (alphas3)==0
						mkmat alpha02, mat(alpha02)
						set seed 280115
						sureg (alpha01) (alpha02), vce(cluster nwbcode)
						test [alpha01]_cons + [alpha02]_cons -1 = 0
						local al_chi2_stat = r(chi2)
						local al_chi2_p = r(p)
						
						*** Run PCDID model ***
						eststo: xtmg y dem_interaction $xvars2 ///
							c1factor*N c2factor*N ///
							if sample==1 & treat_inter==1, res(resid) robust
						mat betas = e(betas)
						estadd scalar treated = e(N_g)
						estadd scalar obs = e(N)
						local treated = e(N_g)
						local obs = e(N)
						mat b1 = e(b)
						mat v1 = e(V)
						scalar b_1 = b1[1,1]
						scalar se_1 = sqrt(v1[1,1])	
						scalar mse = e(sigma)
						estadd scalar mse = e(sigma)
						capture xtcd resid if e(sample)
						capture scalar cdtest = r(pesaran)
						capture estadd scalar cdtest = r(pesaran)
						estadd scalar fac = `f'
						estadd scalar x_test_chi2 = `chi2_stat'
						estadd scalar x_test_p = `chi2_p'
						estadd scalar alpha_b = `al_chi2_stat'
						estadd scalar alpha_p = `al_chi2_p'
						estadd scalar controlN1 = `below'
						estadd scalar controlobs1 = `belowobs'
						estadd scalar controlN2 = `below_too'
						estadd scalar controlobs2 = `belowobs_too'
						*estadd scalar controlN3 = `below_inter'
						*estadd scalar controlobs3 = `belowobs_inter'
						estadd scalar final_year = `kkk'
						estadd scalar start_year = `kk'
						est store model`i'
						
						keep if sample==1 & treat_inter==1
						tabstat nwbcode if sample==1 & treat_inter==1, by(wbcode) save
						tabstatmat nwbcodestat, nototal
						tabstat year if sample==1 & treat_inter==1, by(wbcode) stats(min max) save
						tabstatmat yearstat, nototal
						tabstat dem_dum_count_inter trudem_dum_count_inter dem_dum_exposure_inter ///
							sample_obs_inter autoc_exposure_inter, by(wbcode) stats(mean) save
						tabstatmat exposure, nototal
						drop _all
						
						svmat nwbcodestat, name(nwbcode)
						rename nwbcode1 nwbcode
						svmat yearstat, name(year)
						rename year1 startyear
						rename year2 endyear
						svmat exposure, name(count)
						rename count1 flips
						rename count2 democratisations
						rename count3 exposure
						rename count4 total_obs
						rename count5 nonexposure
						svmat betas, name(beta)
						gen wbcode=""
						replace wbcode="AFG" if nwbcode==	2
						replace wbcode="AGO" if nwbcode==	3
						replace wbcode="ALB" if nwbcode==	5
						replace wbcode="ARE" if nwbcode==	7
						replace wbcode="ARG" if nwbcode==	8
						replace wbcode="ARM" if nwbcode==	9
						replace wbcode="AUS" if nwbcode==	12
						replace wbcode="AUT" if nwbcode==	13
						replace wbcode="AZE" if nwbcode==	14
						replace wbcode="BDI" if nwbcode==	15
						replace wbcode="BEL" if nwbcode==	17
						replace wbcode="BEN" if nwbcode==	18
						replace wbcode="BFA" if nwbcode==	19
						replace wbcode="BGD" if nwbcode==	20
						replace wbcode="BGR" if nwbcode==	21
						replace wbcode="BHR" if nwbcode==	22
						replace wbcode="BIH" if nwbcode==	24
						replace wbcode="BLR" if nwbcode==	25
						replace wbcode="BOL" if nwbcode==	29
						replace wbcode="BRA" if nwbcode==	30
						replace wbcode="BRB" if nwbcode==	31
						replace wbcode="BWA" if nwbcode==	36
						replace wbcode="CAF" if nwbcode==	37
						replace wbcode="CAN" if nwbcode==	38
						replace wbcode="CHE" if nwbcode==	39
						replace wbcode="CHL" if nwbcode==	40
						replace wbcode="CHN" if nwbcode==	41
						replace wbcode="CIV" if nwbcode==	42
						replace wbcode="CMR" if nwbcode==	43
						replace wbcode="COG" if nwbcode==	45
						replace wbcode="COL" if nwbcode==	46
						replace wbcode="COM" if nwbcode==	47
						replace wbcode="CPV" if nwbcode==	48
						replace wbcode="CRI" if nwbcode==	49
						replace wbcode="CUB" if nwbcode==	51
						replace wbcode="CYP" if nwbcode==	52
						replace wbcode="CZE" if nwbcode==	53
						replace wbcode="DEU" if nwbcode==	55
						replace wbcode="DJI" if nwbcode==	56
						replace wbcode="DNK" if nwbcode==	58
						replace wbcode="DOM" if nwbcode==	59
						replace wbcode="DZA" if nwbcode==	60
						replace wbcode="ECU" if nwbcode==	61
						replace wbcode="EGY" if nwbcode==	62
						replace wbcode="ESP" if nwbcode==	64
						replace wbcode="EST" if nwbcode==	65
						replace wbcode="ETH" if nwbcode==	66
						replace wbcode="FIN" if nwbcode==	67
						replace wbcode="FRA" if nwbcode==	69
						replace wbcode="GAB" if nwbcode==	71
						replace wbcode="GBR" if nwbcode==	72
						replace wbcode="GEO" if nwbcode==	73
						replace wbcode="GHA" if nwbcode==	74
						replace wbcode="GIN" if nwbcode==	76
						replace wbcode="GMB" if nwbcode==	77
						replace wbcode="GNB" if nwbcode==	78
						replace wbcode="GNQ" if nwbcode==	79
						replace wbcode="GRC" if nwbcode==	80
						replace wbcode="GTM" if nwbcode==	83
						replace wbcode="HKG" if nwbcode==	87
						replace wbcode="HND" if nwbcode==	89
						replace wbcode="HRV" if nwbcode==	91
						replace wbcode="HTI" if nwbcode==	92
						replace wbcode="HUN" if nwbcode==	93
						replace wbcode="IDN" if nwbcode==	95
						replace wbcode="IND" if nwbcode==	96
						replace wbcode="IRL" if nwbcode==	97
						replace wbcode="IRN" if nwbcode==	98
						replace wbcode="IRQ" if nwbcode==	99
						replace wbcode="ISL" if nwbcode==	100
						replace wbcode="ISR" if nwbcode==	101
						replace wbcode="ITA" if nwbcode==	102
						replace wbcode="JAM" if nwbcode==	103
						replace wbcode="JOR" if nwbcode==	104
						replace wbcode="JPN" if nwbcode==	105
						replace wbcode="KAZ" if nwbcode==	106
						replace wbcode="KEN" if nwbcode==	107
						replace wbcode="KGZ" if nwbcode==	108
						replace wbcode="KHM" if nwbcode==	109
						replace wbcode="KOR" if nwbcode==	112
						replace wbcode="KWT" if nwbcode==	113
						replace wbcode="LAO" if nwbcode==	115
						replace wbcode="LBN" if nwbcode==	116
						replace wbcode="LBR" if nwbcode==	117
						replace wbcode="LBY" if nwbcode==	118
						replace wbcode="LKA" if nwbcode==	120
						replace wbcode="LSO" if nwbcode==	121
						replace wbcode="LTU" if nwbcode==	122
						replace wbcode="LUX" if nwbcode==	123
						replace wbcode="LVA" if nwbcode==	124
						replace wbcode="MAR" if nwbcode==	126
						replace wbcode="MDA" if nwbcode==	128
						replace wbcode="MDG" if nwbcode==	129
						replace wbcode="MEX" if nwbcode==	132
						replace wbcode="MLI" if nwbcode==	135
						replace wbcode="MLT" if nwbcode==	136
						replace wbcode="MMR" if nwbcode==	137
						replace wbcode="MNE" if nwbcode==	138
						replace wbcode="MNG" if nwbcode==	139
						replace wbcode="MOZ" if nwbcode==	140
						replace wbcode="MRT" if nwbcode==	141
						replace wbcode="MUS" if nwbcode==	143
						replace wbcode="MWI" if nwbcode==	144
						replace wbcode="MYS" if nwbcode==	145
						replace wbcode="NAM" if nwbcode==	146
						replace wbcode="NER" if nwbcode==	148
						replace wbcode="NGA" if nwbcode==	149
						replace wbcode="NIC" if nwbcode==	150
						replace wbcode="NLD" if nwbcode==	151
						replace wbcode="NOR" if nwbcode==	152
						replace wbcode="NPL" if nwbcode==	153
						replace wbcode="NZL" if nwbcode==	156
						replace wbcode="OMN" if nwbcode==	158
						replace wbcode="PAK" if nwbcode==	159
						replace wbcode="PAN" if nwbcode==	160
						replace wbcode="PER" if nwbcode==	161
						replace wbcode="PHL" if nwbcode==	162
						replace wbcode="POL" if nwbcode==	165
						replace wbcode="PRK" if nwbcode==	168
						replace wbcode="PRT" if nwbcode==	170
						replace wbcode="PRY" if nwbcode==	171
						replace wbcode="QAT" if nwbcode==	176
						replace wbcode="RUS" if nwbcode==	179
						replace wbcode="RWA" if nwbcode==	180
						replace wbcode="SAU" if nwbcode==	181
						replace wbcode="SDN" if nwbcode==	183
						replace wbcode="SEN" if nwbcode==	184
						replace wbcode="SGP" if nwbcode==	185
						replace wbcode="SLE" if nwbcode==	187
						replace wbcode="SLV" if nwbcode==	188
						replace wbcode="STP" if nwbcode==	195
						replace wbcode="SVK" if nwbcode==	198
						replace wbcode="SVN" if nwbcode==	199
						replace wbcode="SWE" if nwbcode==	200
						replace wbcode="SWZ" if nwbcode==	201
						replace wbcode="SYC" if nwbcode==	203
						replace wbcode="SYR" if nwbcode==	204
						replace wbcode="TCD" if nwbcode==	205
						replace wbcode="TGO" if nwbcode==	206
						replace wbcode="THA" if nwbcode==	207
						replace wbcode="TJK" if nwbcode==	208
						replace wbcode="TKM" if nwbcode==	209
						replace wbcode="TTO" if nwbcode==	214
						replace wbcode="TUN" if nwbcode==	215
						replace wbcode="TUR" if nwbcode==	216
						replace wbcode="TZA" if nwbcode==	220
						replace wbcode="UGA" if nwbcode==	221
						replace wbcode="UKR" if nwbcode==	222
						replace wbcode="URY" if nwbcode==	223
						replace wbcode="USA" if nwbcode==	224
						replace wbcode="UZB" if nwbcode==	225
						replace wbcode="VEN" if nwbcode==	228
						replace wbcode="VNM" if nwbcode==	229
						replace wbcode="YEM" if nwbcode==	235
						replace wbcode="ZAF" if nwbcode==	238
						replace wbcode="ZMB" if nwbcode==	240
						replace wbcode="ZWE" if nwbcode==	241
						order nwbcode wbcode startyear endyear beta1 flips democratisations exposure total_obs nonexposure
						keep nwbcode wbcode startyear endyear beta1 flips democratisations exposure total_obs nonexposure
						gen regime_def = "inter"
						local startyearx "$startyear"
						local endyearx "$endyear" 
						local democracyx "$democracy" 
						local democracyxx "$democracy_too" 
						save "$outputfolder/PCDID_Maddison2021_Drilling_st_Inter_`startyearx'_`endyearx'_`democracyx'_`democracyxx'_`f'.dta", replace
						
						local i = `i'+1

					}
					* end of quietly
												
					di in ye _n "*****************************************************************"
					di in ye "*** Single regime cutoff defined by " in gr "dem_$democracy" 
					di in ye "*****************************************************************"
					di in ye "*** PCDID (MG) estimate with" in gr " `f' " in ye "common factor(s) included      ***"
					di in ye "*****************************************************************"
					di in ye "Nominal sample period: " in gr "$startyear" in ye " to " in gr "$endyear "
					di in ye "resulting in " in gr "`treated'" in ye " treated countries with " ///
						in gr "`obs'" in ye " observations. "
					di in ye "MSE: " in gr %5.3f mse
					di in ye "Control sample 1: " in gr "`below'" in ye " countries with " in gr "`belowobs'" in ye " obs."
					di in ye "Control sample 2: " in gr "`below_too'" in ye " countries with " in gr "`belowobs_too'" in ye " obs."
					*di in ye "Control sample: " in gr "`below_inter'" in ye " countries with " in gr "`belowobs_inter'" in ye " obs."
					di in ye "*****************************************************************"
					*di in ye "Estimates for 1st democracy indicator: " in gr "dem_$democracy"
					*di in ye "Rob. Beta: " in gr %5.3f b_1 in ye "  SE: " in gr %5.3f se_1 ///
					*	in ye "  t-ratio: " in gr %5.3f b_1/se_1 in ye "  p: " in gr %5.3f 2*(normal(-abs(b_1/se_1)))
					*di in ye "*****************************************************************"
					di in ye "Estimates for interaction with: "  in gr "dem_$democracy_too"
					di in ye "Rob. Beta: " in gr %5.3f b_1 in ye "  SE: " in gr %5.3f se_1 ///
						in ye "  t-ratio: " in gr %5.3f b_1/se_1 in ye "  p: " in gr %5.3f 2*(normal(-abs(b_1/se_1)))
					di in ye "*****************************************************************"
					di in ye "*** Alpha Test: " in gr %5.2f `al_chi2_stat' in ye " (" in gr %5.3f `al_chi2_p' in ye ").      ***"
					di in ye "*****************************************************************"
					di in ye "*** Test of controls " in gr %5.2f `chi2_stat' in ye " (" in gr %5.3f `chi2_p' in ye ").      ***"
					di in ye "*****************************************************************"
					local ii = `i'-1
					di in red "`ii'"
					
					restore	
					
					}
					* end of the factor loop
				}
				* end of second dem indicator loop
			}
			* end of first dem indicator loop 
		}
		* end of start year loop
	}
	* end of end year loop

	local mytitle "PCDID"
	local esttabformat  "nodepvars lines compress label nogaps numbers star(* 0.10 ** 0.05 *** 0.01)"
	local esttabstats "stats(start_year final_year fac treated obs controlN1 controlobs1 controlN2 controlobs2 x_test_chi2 x_test_p alpha_b alpha_p, fmt(0 0 0 0 0 0 0 2 2 3 3))"
	local esttabkeep "keep($demo_list2 $demo_list_2 dem_interaction $xvars2) order($demo_list2 $demo_list_2 dem_interaction $xvars2)"
	esttab model* using "$outputfolder/20240717_Maddison2021_PCDID_Drilling_`democracyx'_`democracyxx'_alt_model_all_factors.csv", b(3) se(3) abs br ///
		`esttabformat' `esttabkeep' `esttabstats' style(tab) replace nodep

}
* end of the loop related to the variation of the threshold (-.25sd -.125sd 0 [mean] +.125sd +.25sd)

exit

***
*** Graph: Comparison interaction with poly and lib 
***

use "$outputfolder/PCDID_Maddison2024_Drilling_1959_2018_st_multiple_poly_0_diagnostics_4.dta", clear 
rename beta1 beta
drop regime_def
local z "startyear endyear beta flips democratisations exposure total_obs nonexposure"
foreach k of local z{
	rename `k' `k'_poly
}
drop nwbcode
merge wbcode using "$outputfolder/PCDID_Maddison2024_Drilling_1959_2018_st_multiple_lib_0_diagnostics_4.dta", sort 
rename beta1 beta
drop regime_def
local z "startyear endyear beta flips democratisations exposure total_obs nonexposure"
foreach k of local z{
	rename `k' `k'_lib
}
drop nwbcode _merge
merge wbcode using "$outputfolder/PCDID_Maddison2021_Drilling_st_Inter_1959_2018_lib_0_poly_0_4.dta", sort
rename beta1 beta
drop regime_def
local z "startyear endyear beta flips democratisations exposure total_obs nonexposure"
foreach k of local z{
	rename `k' `k'_libpoly
}

tab startyear_poly, gen(by_)
tab startyear_lib, gen(b_y_)
tab startyear_libpoly, gen(be_y__)

*** No horserace
mrunning beta_poly exposure_poly flips_poly if !missing(beta_poly,beta_lib, beta_libpoly), gen(change_) ci gense(changese_)  adjust(by_*) nograph
rename change_1 change_poly
rename changese_1 changese_poly
order change_* changese_*
drop change_2  changese_2
drop  change_3-change_18  changese_3-changese_18
gen tstat_poly = change_poly/changese_poly

mrunning beta_lib exposure_lib flips_lib if !missing(beta_poly,beta_lib, beta_libpoly), gen(change_) ci gense(changese_) adjust(b_y_*) nograph
rename change_1 change_lib
rename changese_1 changese_lib
order change_* changese_*
drop change_2  changese_2
drop change_3-change_16  changese_3-changese_16
gen tstat_lib = change_lib/changese_lib

mrunning beta_libpoly exposure_libpoly flips_libpoly if !missing(beta_poly,beta_lib, beta_libpoly), gen(change_) ci gense(changese_)  adjust(by_*) nograph
rename change_1 change_libpoly
rename changese_1 changese_libpoly
order change_* changese_*
drop change_2  changese_2
drop  change_3-change_18  changese_3-changese_18
gen tstat_libpoly = change_libpoly/changese_libpoly

rreg change_libpoly if !missing(change_lib,change_poly)
rreg change_lib  if !missing(change_libpoly,change_poly)
rreg change_poly  if !missing(change_lib,change_libpoly)

twoway ///
	(hist exposure_lib, yaxis(2) w(5) start(0) freq fcolor(gs9%50) lcolor(white) barwidth(4.5)) ///
	(line change_poly exposure_poly, sort ///
		lw(.8) lcolor(pink*1.5) lpattern(shortdash)) ///
	(scatter change_poly exposure_poly if abs(tstat_poly)<1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(medlarge) mfcolor(white%60) mlcolor(pink*1.5%60) mlwidth(thin)) ///
	(scatter change_poly exposure_poly if abs(tstat_poly)>1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(medlarge) mfcolor(pink%60) mlcolor(pink*1.5%60) mlwidth(thin)) ///
	(line change_lib exposure_lib, sort ///
		lw(.8) lcolor(midblue*1.5) lpattern(dash)) ///
	(scatter change_lib exposure_lib if abs(tstat_lib)<1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(medlarge) mfcolor(white%60) mlcolor(midblue*1.5%60) mlwidth(thin)) ///
	(scatter change_lib exposure_lib if abs(tstat_lib)>1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(medlarge) mfcolor(midblue%60) mlcolor(midblue*1.5%60) mlwidth(thin)) ///
	(line change_libpoly exposure_libpoly, sort ///
		lw(.8) lcolor(emerald*1) lpattern(solid)) ///
	(scatter change_libpoly exposure_libpoly if abs(tstat_libpoly)<1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(large) mfcolor(white%60) mlcolor(emerald*1.25%60) mlwidth(thin)) ///
	(scatter change_libpoly exposure_libpoly if abs(tstat_libpoly)>1.645, sort jitterseed(280115) jitter(2) ///
		 msymbol(O)  msize(large) mfcolor(emerald%60) mlcolor(emerald*1.25%60) mlwidth(thin)) ///
, xsize(9) ysize(4) scheme(s2mono) yline(0, lpattern(shortdash) lw(.4) lcolor(black)) ///
	xtitle("Years Spent in Democracy (Length of Treatment)", size(medsmall)) ///
	ytitle("Democracy Effect (in %)", size(medsmall))  ///
	ytitle("Country Count", size(medsmall) axis(2))  ///
	yscale(range(-10.1 40.1))   ///
	ylabel(-10(10)40, labsize(small) glcolor(gs8) grid glwidth(.2) glp(shortdash)  angle(0)) ///
	ylabel(0(2)12, labsize(small) glcolor(gs8) angle(0) axis(2)) ///
	xlabel(0(5)60, labsize(small)) graphregion(color(white)) ///
	xmtick(##5) ymtick(##2, axis(2))  ///
	legend(rows(2) order(- "Democracy Effect:" 2 5 8 - "Countries:" 1 - "Significant (10% level):" 4 7 10) ///
	label(10 "") label(2 "Polyarchy > mean") label(5 "Liberal Component > mean") label(8 "Interaction")	 ///
	label(4 "")	label(7 "") label(1 "64")	size(medsmall)  symxsize(*.5) holes(11)) 
*graph export "$texfolder/PCDID_Maddison2021_mid_levels_st_interaction_new.png", as(png) replace width(1200)
graph export "$texfolder/PCDID_Maddison2024_mid_levels_st_interaction_new.eps", as(eps) replace


***
*** Graph: Comparison interaction_0 with other cutoffs 
***

use "$outputfolder/PCDID_Maddison2021_Drilling_st_Inter_1959_2018_lib_neg25_poly_neg25_4.dta", clear
rename beta1 beta
drop regime_def
local z "startyear endyear beta flips democratisations exposure total_obs nonexposure"
foreach k of local z{
	rename `k' `k'_neg25
}
drop nwbcode
local zz "neg125 0 125 25"
foreach kk of local zz{	
	merge wbcode using "$outputfolder/PCDID_Maddison2021_Drilling_st_Inter_1959_2018_lib_`kk'_poly_`kk'_4.dta", sort 
	rename beta1 beta
	drop regime_def _merge
	local z "startyear endyear beta flips democratisations exposure total_obs nonexposure"
	foreach k of local z{
		rename `k' `k'_`kk'
	}
}

local z "neg25 neg125 0 125 25"
foreach k of local z{
	tab startyear_`k', gen(by_) 
	mrunning beta_`k' exposure_`k' flips_`k', gen(change__) ci gense(changese__)  adjust(by_*) nograph
	rename change__1 change_`k'
	rename changese__1 changese_`k'
	drop change__* changese__*
	gen tstat_`k' = change_`k'/changese_`k'
	drop by_*
}

twoway ///
	(line change_neg25 exposure_neg25, sort /// 1
		lcolor(gs2%60) lwidth(.8) lpattern(solid)) ///
	(scatter change_neg25 exposure_neg25 if abs(tstat_neg25)<1.645, sort jitterseed(280115) jitter(2) ///2
		 msymbol(O)  msize(medlarge) mfcolor(white%40) mlcolor(gs2%60) mlwidth(thin)) ///
	(scatter change_neg25 exposure_neg25 if abs(tstat_neg25)>1.645, sort jitterseed(280115) jitter(2) ///3
		 msymbol(O)  msize(medlarge) mfcolor(gs2%60) mlcolor(gs2%60) mlwidth(thin)) ///
	(line change_neg125 exposure_neg125, sort ///4
		lcolor(gs5%60) lwidth(.8) lpattern(solid)) ///
	(scatter change_neg125 exposure_neg125 if abs(tstat_neg125)<1.645, sort jitterseed(280115) jitter(2) ///5
		 msymbol(O)  msize(medlarge) mfcolor(white%40) mlcolor(gs5%60) mlwidth(thin)) ///
	(scatter change_neg125 exposure_neg125 if abs(tstat_neg125)>1.645, sort jitterseed(280115) jitter(2) ///6
		 msymbol(O)  msize(medlarge) mfcolor(gs5%60) mlcolor(gs5%60) mlwidth(thin)) ///
	(line change_125 exposure_125, sort ///7
		lcolor(gs8%60) lwidth(.8) lpattern(solid)) ///
	(scatter change_125 exposure_125 if abs(tstat_125)<1.645, sort jitterseed(280115) jitter(2) ///8
		 msymbol(O)  msize(medlarge) mfcolor(white%40) mlcolor(gs8%60) mlwidth(thin)) ///
	(scatter change_125 exposure_125 if abs(tstat_125)>1.645, sort jitterseed(280115) jitter(2) ///9
		 msymbol(O)  msize(medlarge) mfcolor(gs8%60) mlcolor(gs8%60) mlwidth(thin)) ///
	(line change_25 exposure_25, sort ///10
		lcolor(gs11%60) lwidth(.8) lpattern(solid)) ///
	(scatter change_25 exposure_25 if abs(tstat_25)<1.645, sort jitterseed(280115) jitter(2) ///11
		 msymbol(O)  msize(medlarge) mfcolor(white%40) mlcolor(gs11%60) mlwidth(thin)) ///
	(scatter change_25 exposure_25 if abs(tstat_25)>1.645, sort jitterseed(280115) jitter(2) ///12
		 msymbol(O)  msize(medlarge) mfcolor(gs11%60) mlcolor(gs11%60) mlwidth(thin)) ///
	(line change_0 exposure_0, sort yaxis(2) ///13
		lcolor(emerald*1.25) lwidth(.8) lpattern(solid)) ///
	(scatter change_0 exposure_0 if abs(tstat_0)<1.645, sort jitterseed(280115) jitter(2) ///14
		 msymbol(O)  msize(large) mfcolor(white%40) mlcolor(emerald*1.5) mlwidth(thin)  yaxis(2)) ///
	(scatter change_0 exposure_0 if abs(tstat_0)>1.645, sort jitterseed(280115) jitter(2) ///15
		 msymbol(O)  msize(large) mfcolor(emerald*1) mlcolor(emerald*1.5) mlwidth(thin)  yaxis(2)) ///
, xsize(9) ysize(4) scheme(s2mono) yline(0, lpattern(shortdash) lw(.4) lcolor(black)) ///
	xtitle("Years Spent in Regime (Length of Treatment)", size(medsmall)) ///
	ytitle("Democracy Effect (in %)", size(medsmall))  ///
	yscale(range(-10.1 62.1))   ///
	ylabel(-10(10)60, labsize(small) glcolor(gs8) angle(0) axis(2)) ///
	ylabel(-10(10)60, labsize(small) glcolor(gs8) grid glwidth(.2) glp(shortdash)  angle(0)) ///
	xlabel(0(5)60, labsize(small)) graphregion(color(white)) ///
	ytitle("Democracy Effect (in %)", size(medsmall) axis(2))  ///
	yscale(range(-10.1 62.1) axis(2)) 	xmtick(##5) ///
	legend(rows(2) order(- "Interaction for Cut-offs at" 1 4 13 7 10 - "Significant (10% level)" ///
	3 6 15 9 12) ///
	label(1 "-1/4 sd")  label(4 "-1/8 sd") label(13 "mean") label(7 "+1/8 sd") label(10 "+1/4 sd") ///
	label(3 "") label(6 "") label(15 "") label(9 "") label(12 "") ///
	size(medsmall) symxsize(*.5)) 
*graph export "$texfolder/PCDID_Maddison2021_mid_levels_stvarious_interaction.png", as(png) replace width(1200)
graph export "$texfolder/PCDID_Maddison2024_mid_levels_stvarious_interaction.eps", as(eps) replace



*** Regression results (comparison)
cls

estimates clear 
local i = 1

local zz "neg25 neg125 0 125 25"
foreach kk of local zz{
	quietly{
		use "$outputfolder/PCDID_Maddison2024_Drilling_1959_2018_st_multiple_poly_`kk'_diagnostics_4.dta", clear 
		rename beta1 beta
		drop regime_def
		local z "startyear endyear beta flips democratisations exposure total_obs nonexposure"
		foreach k of local z{
			rename `k' `k'_poly
		}
		drop nwbcode
		merge wbcode using "$outputfolder/PCDID_Maddison2024_Drilling_1959_2018_st_multiple_lib_`kk'_diagnostics_4.dta", sort 
		rename beta1 beta
		drop regime_def
		local z "startyear endyear beta flips democratisations exposure total_obs nonexposure"
		foreach k of local z{
			rename `k' `k'_lib
		}
		drop nwbcode _merge
		merge wbcode using "$outputfolder/PCDID_Maddison2021_Drilling_st_Inter_1959_2018_lib_`kk'_poly_`kk'_4.dta", sort
		rename beta1 beta
		drop regime_def
		local z "startyear endyear beta flips democratisations exposure total_obs nonexposure"
		foreach k of local z{
			rename `k' `k'_libpoly
		}
	}
	di in gr "**********************************************************************************"
	di in ye "Estimates for `kk'"
	di in gr "**********************************************************************************"


//	eststo: rreg beta_poly if !missing(beta_poly, beta_lib, beta_libpoly)	
//	eststo: rreg beta_lib if !missing(beta_poly, beta_lib, beta_libpoly)	
 	eststo: rreg beta_libpoly if !missing(beta_poly, beta_lib, beta_libpoly)

	est store model_`i'

	
// 	eststo: reg total_obs_poly if !missing(beta_poly, beta_lib, beta_libpoly)	
//	eststo: reg total_obs_lib if !missing(beta_poly, beta_lib, beta_libpoly)
	eststo: reg total_obs_libpoly if !missing(beta_poly, beta_lib, beta_libpoly)

	est store modl_`i'

// 	eststo: qreg exposure_poly if !missing(beta_poly, beta_lib, beta_libpoly)
// 	eststo: qreg exposure_lib if !missing(beta_poly, beta_lib, beta_libpoly)
 	eststo: qreg exposure_libpoly if !missing(beta_poly, beta_lib, beta_libpoly)

	est store mod_`i'

	local i = `i' + 1	


}


local esttabformat  "lines label nogaps numbers star(* 0.10 ** 0.05 *** 0.01) staraux mlabels(, depvars) note("") nolegend"
	
esttab model_* , b(3) se(3) abs `esttabformat' br

esttab modl_* , b(3) abs `esttabformat'

esttab mod_* , b(0) abs `esttabformat'
